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We consider here the recently proposed closed form formula in terms of the Meijer G-functions 
for the probability density functions g a (x) of one-sided bevy stable distributions with rational index 
a — l/k, with < a < 1. Since one-sided bevy and Mittag-befner distributions are known to be 
related, this formula could also be useful for calculating the probability density functions p a {x) of 
the latter. We show, however, that the formula is computationally inviable for fractions with large 
denominators, being unpractical even for some modest values of I and k. We present a fast and 
accurate numerical scheme, based on an early integral representation due to Mikusinski, for the 
evaluation of g a (x) and pa(x), their cumulative distribution function and their derivatives for any 
real index a £ (0, 1). As an application, we explore some properties of these probability density 
functions. In particular, we determine the location and value of their maxima as functions of the 
index a. We show that a ~ 0.567 and a ~ 0.605 correspond, respectively, to the one-sided bevy 
and Mittag-befner distributions with shortest maxima. We close by discussing how our results can 
elucidate some recently described dynamical behavior of intermittent systems. 

PACS numbers: 05.40.Fb, 02.50.Ng, 02.60.Jh 



I. INTRODUCTION 

One-sided Levy stable distributions [l|, [|| are ubiqui- 
tous in many modern research areas where quantitative 
and statistical analysis play a major role. (For recent 
reviews, see, besides 0, the references of [||.) The prob- 
ability density function of one-sided Levy distribution of 
index a, g a (x), can be defined by means of its Laplace 
transform as [l| 



e bX g a (x) dx = exp (— s a ) . 



(1) 



for s > 0, with < a < 1. Unfortunately, in spite 
of its broad applicability, exact solutions of Eq. (JT|) are 
available only for a few particular values of a. (See, for 
instance, the Appendix A of [![. We notice also that 
there are some available Mathematica [| and Matlab [| 
packages for the numerical evaluation of g a (x).)ln this 
context, the recent work of Penson and Gorska Q is cer- 
tainly interesting and relevant since they describe a for- 
mal solution of Eq. (JTJ) for any rational a. In fact, they 
show that a formula presented without proof in a table 
of inverse Laplace transforms Q could be used to write 
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A(Z,0) 
A(fc,0) 
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where G™^" (z ^ ^ is the Meijer G-function 
A(fc, a) is the list of k elements given by 

a a + 1 a + k — 1 
fc' k k 



A(k,a) 



(2) 
and 

(3) 



We consider the formula ([2]) an important advance. Since 
the Meijer G-function is available in several computer al- 
gebra systems, the function gi/k{x) could in principle be 
evaluated with little programming effort. We notice that 
the restriction to rational values of a in Eq. ([2} does not 
represent any real problem here. As we will see below, 
the function g a (x) is continuous in a and, hence, one 
might compute from Eq. ([2]) a rational a approxima- 
tion for (x) with any prescribed accuracy. Penson and 
Gorska |3| use Eq. ^ to derive other series expression 
for gi/k{x) and to infer some of its properties. Certainly, 
the mathematical literature about the Meijer G-function 
(see, for instance, [8| and the references therein) will be 
extremely valuable for the derivation of many other prop- 
erties of gi/k{x) denned by Eq. (|2|). 

Furthermore, since one-sided Levy and Mittag-Lcfflcr 
distributions are known to be related the formula 
([2]) is also relevant for calculating the probability density 

functions pa^ (x) of Mittag-Leffler distributions with ra- 
tional index a. We recall that pa\x) is also defined from 
its Laplace transform as well, 



e- sx p^\x)dx 
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(4) 



for s > 0, with < a < 1. The right-handed side of 
Eq. §4$ is a particular case of the so-called Mittag-Leffler 
function [8(, which reduces to the usual exponential for 
a = 1. The free parameter r can be fixed, for instance, 
by demanding a given first moment for pa\x). In par- 
ticular, since we have from Eq. (jH that 



p^(qx) = q- 1 p^ 



(5) 
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for any q > 0, one can assume r = 1 without loss of gen- 
erality. In this case, the superscript is simply dropped. 
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The respective cumulative distribution functions associ- 
ated to p a (x) and g a (x) are known to be related by |9[ 



e Q (x) = i-A Q (x- 1 /") 



which leads to 



(6) 



(7) 



The relation ([7]) allows the computation of p a {x) by 
means of the Meijer G-function for rational a, thanks 
to the Penson and Gorska formula ([2]). This is a consid- 
erable advance since, as in the previous case, no closed 
form solution of Eq. (j4]) is known. 

However, the condensed and apparently simple form 
of Eq. @ hides a practical pitfall. The evaluation of 
Eq. ([2]) is computationally viable only for modest values 
of k and For instance, by using the Maple procedure 
provided by Penson and Gorska [lo| . we can plot the 
graphics of {72/3 (x) for x € [0,2] instantaneously in an 
Intel Core i7 computer running Maple version 14. In 
order to generate the same graphics with, for instance, 
l/k = 20/31, some CPU minutes are necessary. For 
l/k = 200/301, we need almost a half an hour to evaluate 
a single value of gi/k(x)\ We could not evaluate Eq. ([2]) 
for l/k = 2000/3001 in any reasonable amount of time. 
Mathematica presents a similar performance. These re- 
strictions, obviously, jeopardize the practical utility of 
expression ([2]) since one cannot calculate in reasonable 
time good approximations to the one-sided Levy distri- 
bution for any a. One can understand the rapidity with 
the evaluation of Eq. @ becomes unpractical when the 
values of I and k increase by recalling the definition of 
the Meijer G-function H 



rim , n 
S,9 
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s ds, 



where L is a carefully chosen integration path on the 
complex plane. It is possible also to write the Meijer 
G-function as a sum of m terms involving T function 
products as those ones of the integrand in Eq. flSJ) and 
the generalized hypergeometric functions P F q —i [8|. As 
one can see, when asking Maple to evaluate Eq. @ for 
l/k = 2000/3001, one is basically demanding the evalua- 
tion of an integral with more than five thousands T func- 
tion terms in the integrand, or an intricate combination 
of more than three thousands generalized hypergeometric 
functions! Hence, it is not a surprise to have a consider- 
able performance degradation for large values of I and k. 
Another problem with the Maple procedure based in Eq. 
^ is that it does not deal efficiently with reducible frac- 
tions. For instance, Maple is not able to reduce 35/10 (x) 
to <7i/2(x). Moreover, the numerical evaluation of the 
former is much more time and memory consuming than 
the latter. 



The purpose of the present work is to show that one 
can compute numerically, in an effective and efficient way, 
the probability density functions g a (x) and p a (x) with 
arbitrary real index a € (0,1). Our start point is the 
Mikusinski's integral representation for g a {x) (Tl| 



la{x) 



a 1 
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with < a < 1, where 



ue u dip, 
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(9) 



(10) 



The integral representation (j9|) has already proven its 
relevance. In fact, Mikusinski used it to derive more than 
40 years ago the some very useful asymptotic expressions 
for g a (x), namely 



g a (x) « K 



exp(-Ax~ Q /( 1 - Q )) 



(11) 



(12) 



2 ,(2-Q)/(2-2a) ' 

valid for x — > + and 

5ct (x) « Mx-( 1+Q \ 
valid for x — > 00, where 

A = (l-a)a! a/(1 ~ a) , (13) 

l/(2-2a) 

(14) 



K = 



v^a-a) 



smmr . . 
M = r(l + a). 



(15) 



It is easy to check from Eq. (|9]) that g a (x) is non- negative 
and smooth in a 6 (0, 1). Due to Eq. (0, and one has 
the following asymptotic behavior for p a (x) 



Pa{x) -> r(l + a), 

air 



(16) 



for x — > + and 



pa ( x ) „ ^ a; (2«-l)/(2-2 Q ) exp (- Ax im-<x^ ( (17) 

valid for x — > 00. In the next section, we show how to 
use Eq. ^ to evaluate numerically one-sided Levy and 
Mittag-Leffler probability densities, their cumulative dis- 
tribution function, and their derivatives in a very efficient 
and reliable way. 

We notice that one-sided Levy stable distributions can 
be alternatively expressed by means of Fox H-functions 
, which are a further generalization of the Meijer G- 
functions © , and also by means of Wright functions [l3[ . 
Unfortunately, the current knowledge about the analyt- 
ical structure of these functions is still little developed. 
We wish also to stress here that the numerical compu- 
tation of stable distributions is not a new problem and 
several algorithms are already available in the literature 
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and even commercially. In particular, Nolan pro- 
posed a robust algorithm based on the integration of the 
so-called Zolotarev's (M) representation for stable distri- 
butions, which is the base of Mathematica [f| and Matlab 
@ packages. An updated reference list on the subject 
can be found in (la ]. However, as we will see, Mikusin- 
ski's representation ^ allows the numerical evaluation 
of one-sided Levy stable and Mittag-LefHer distributions 
with little programming and computational efforts and 
with the same accuracy of these specialized packages. 
Furthermore, from the Mikusinski's representation one 
will be able to derive some asymptotic expressions with 
special relevance to physical applications. 



II. THE ALGORITHM 

The Mikusinski's integral representation ([9]) involves a 
simple proper integral of a smooth function on the in- 
terval tp € [0, 7r]. Furthermore, it is easy to obtain from 
Eq. (0) some analogous formulas for the derivatives of 
g a (x) and its cumulative distribution function A a (x). In 
particular, we have 



A a (x) 

9a{ x ) 

and 



- / e-^dip, 

7T ' 



1 - a 
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(18) 

1 9a(x) 
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(20) 



where u is given by Eq. ()10j) . The formulas for the 
Mittag-LefHer case can be obtained directly from Eqs. 
© and ([7]). Fig. 1 depicts the integrand in Eq. © for 
some typical values of a and x. The integrands for the 
derivatives of g a (x) and for the Mittag-LefHer case have 
similar aspects. They are all of the type /„ = u n e~ u on 
the interval ip € [0,tt], where u is given by Eq. (fTO)) . It 
is easy to show that: /„ = u[}e _tt(l and f n = for (p = 0, 
with u = (1 - ^(a/x)"^ 1 -^; /„ -)• for p tt; and 
that /„ has a maximum for ip such that u = n, pro- 
vided uq < n. Although the position of such maximum 
does depend on a and x, its value (/„ = n n e~ n ) depends 
only on n. For a given a, larger values of x displace the 
maximum towards p = n, while smaller values does to- 
wards tp — 0. If Uq > n, the unique maximum of /„ is 
at tp = 0. For the cumulative distribution function (fT8"|). 
the integrand corresponds to n = 0. In particular its 
maximum is located at tp — 0, with /q = 1, irrespective 
of the values of x > and < a < 1 . All these functions 
are well-behaved on the interval p € [O^ 71 "] and, conse- 
quently, integrals like Eqs. ©, (TTO| . and can 
be evaluated numerically without major problems. 
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FIG. 1: The integrand / = ue~" of the Mikusinski's represen- 
tation for some values of a and x in the interval tp £ [0, 7r] . 
The curves (a)-(d) correspond, respectively, to the following 
values of (a, x): (0.6,1.0), (0.2,0.1), (0.5,0.2), and (0.5,1.5). 



We have set up an adaptive integration scheme based 
on the publicly available DQAGS routine of SLATEC 
We could integrate Eqs. ©, (H]), (HU), and (JUl) with 
little computational effort demanding a relative error in 
DQAGS smaller than 10 -8 , which is typically attained with 
about 10 iterations of the global adaptative scheme of 
the routine. Our FORTRAN code, available at [l7j], has 
demonstrated to be extremely robust and reliable. In 
order to test it, we have used the case corresponding to 
the Levy distribution with a — 1/2, for which an ex- 
plicit form for the probability density function is known, 
namely the so-called Smirnov's distribution 



S{x) = gi/2(x) 



-1 /Ax 



2y/^X 3 / 2 ' 



(21) 



Fig. [2] shows the functions gi/2(x), Ki/ 2 {x), g' 1 ^ 2 (x), and 
g'^ 2 (x) evaluated numerically with our code. As we see, 
we can calculate g a {x) with very good accuracy and in 
an efficient way. The corresponding data (500 points for 
each curve) for plots like those ones depicted in Fig. [2] are 
generated instantaneously in a Intel Core i7 computer. 
The relation ([7]) and the Smirnov's distribution (|2"Tj) allow 
us to test also the Mittag-LefHer case since they imply 
that 



Pl/2{x) 



(22) 



Our numerical procedure works with similar accuracy for 
this particular Mittag-LefHer probability density, with 
the usual caveats related to extremely small values of x in 
(|22|) . which correspond to large values of x in 31/2(2;) ac- 
cording to (UJ). We also checked the good accuracy of our 
algorith m by comparing the output with Nolan's STABLE 
package [la] . 

The numerical evaluation of the probability densities 
for extreme values of x and a is quite delicate due to 
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FIG. 2: Plots for g 1/2 (x) (a), 2A 1/2 (x) (b), (l/15)gi /a (ar) 
(c), and (l/500)g" /2 (a;) (d), for < x < 1, calculated by 
the numerical integration of Eqs. (J9]), (fT8)> . (fl9|) . and ff20|) 
by means of the SLATEC [3] adaptative integration routine 
DQAGS. The relative errors for all the four curves, calculated 
with respect to the exact Smirnov's distribution (|21l) in the 
interval < x < 5, are smaller than 3 x 10~ 8 . (Numerical 
code available at [171].') 
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FIG. 3: Location and value of the maxima of the probability 
density functions of Levy (curves (a) and (b)) and Mittag- 
Leffler (curves (c) and (d)) distributions. Curves (a) and (c) 
correspond to the location x,(a) of the maxima (left scale), 
as function of a £ (0,1) (horizontal axis). Curves (b) and 
(d) are the value of the maxima (right scale) as function of 
a. The marked points correspond to the shortest maxima 
for each distribution, (a = 0.567, g a (x,) = 0.888) and (a = 
0.605, p a (x„) = 0.509). 



convergence and roundoff problems. For a fixed a and 
x — > and x — > oo, the asymptotic formulas (fTTT). (fl~2|) . 
(dnj), and (TTT1) can be indeed used to estimate the proba- 
bility densities. For fixed x and a very close to and 1, 
other asymptotic expressions are necessary. For small a, 
we have from Eq. (fTTjf 

U« l + ama+ ( In — % —) a. (23) 

\ x sin ip tan p J 

For a fixed x £ (a, 1/a), we have u ~ 1 + a In a for 
small enough a and < p < n, leading to ue~ u ~ 
er l (l — (alna) 2 /2) for small a and < ip < ir. Ap- 
plying these results in Eq. ©, one has 

g a (x) ps — , (24) 

ex 

valid for small a and a < x < 1/a. For x < a, the 
approximation is still valid for small a. In particular, 
we always have g a (x) — ¥ for x — ¥ 0, irrespective of the 
value for a. Since the hypothesis of a < x < 1/a was 
explicitly used, the approximation (f24|) is not supposed 
to be accurate for x — > oo. In this case, Eq. (fl"2")) is the 
correct asymptotic expression for g a {x). For a close to 
1, the situation is a little bit more involved. Introducing 
1 — a = e > 0, we have from Eq. (fT0|) 

a; / s sin 99 

valid for e « 0, < a; < 00, and < y> < 7r. For x > 1, 
u — > for small e, implying that g\- e {x) — > 0. For x < 1, 
u — > 00, also implying gri_ e (x) — > 0. Since, according 
to Eq. |T]), g a (x) is supposed to be normalized for any 



value of a, g a {x) —¥ 00 for x — > 1 and a 1. Hence, 
for a close to 1, g a (x) should be strongly peaked around 
x = 1, resembling an approximation for a ^-function. 
Such behavior could also be inferred by considering the 
limit a — > 1 directly in Eq. ([!]). 



A. The maxima of the distributions 

As an application of our numerical procedures, we will 
explore some properties of the distribution g a {x) and 
p a {x). The location of the maxima of these probabil- 
ity density functions is certainly pertinent to the under- 
standing of the statistical processes governed by them. 
Let us consider fist the case of <? Q (x). The condition 
determining the location x* (a) of the maximum of the 
probability density is, of course, g' a (x*) = 0. From the ap- 
proximations discussed in the preceding section, we have 
that < x*(a) < 1 for < a < 1 and that g a {x*) — > 00 
for a — > and for a — > 1. The zero of g' a (x) can be 
localized in the interval (0, 1) with a prescribed accuracy 
by using, for instance, a simple bisection method. Since 
we have a procedure to calculate 3q(x), one could even 
implement a refinement for the determination of x*(a) 
based, for instance, in Newton-Rapson method. Fig. [3] 
shows the values of x*(a) and g a (x*) for < a < 1. 
Notice that, as expected, we have that g a {x*) — > 00 for 
a — > and a — > 1, in agreement with the approxima- 
tions of last section. The minimal value of g a {x^) is at- 
tained when a = 0.567, corresponding to the one-sided 
a-stable Levy distribution with shortest maximum, for 
which 5 Q (x„) w 0.888. (Fig. gj). 
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FIG. 4: Levy and Mittag-Leffler probability densities with 
shortest maxima (respectively, curves (a), with a = 0.567, 
and (c), with a = 0.605, see Fig. [3J) and their respective 
cumulative distribution function (curves (b) and (d)) in the 
interval < x < 3. 

The situation for p a {x) is rather more involved. For 
a « 1, p a {x) is similar to g a (x), both resembling approx- 
imations of a 8- function around x = 1. Such behavior 
for the Mittag-Leffler case can also be inferred directly 
from the definition (j?]), by considering the limit a — > 1. 
However, in contrast with the previous case, for a ~ 0, 
Pa(x) ~ e~ x , as one can also see by evaluating the limit 
a — > in Eq. Hence, for small a, the maximum of 
p a (x) located at x = and is given by p a (0) ~ 1. In 
fact, we could verify numerically that for a < 1/2, the 
maximum of p a (x) is is always at x — and is given 
by Eq. ([TrJ|. For a > 1/2, we have p' Q (0) > and the 
function attains a maximum for x > and then 

decays. Curiously, as a increases, x*(a) also increases 
and even exceed x = 1, and then return to i = 1, but 
from the right-handed side. This behavior, which will be 
crucial for the discussion of the next section, is depicted 
by the curve (c) in Fig. [3] The Levy and Mittag-Lcfficr 
distribution with shortest maxima are plotted in Fig. 01 



III. DISTRIBUTION OF LYAPUNOV 
EXPONENTS IN INTERMITTENT SYSTEMS 

We can also apply our numerical procedures to eluci- 
date some dynamical problems of physical interest. This 
is the case, for instance, of the distribution of Lyapunov 
exponents in intermittent systems such as the Pomeau- 
Manneville maps Xt+i = x t + ax\ (mod 1) considered re- 
cently in [l9j]. For z > 2, theses systems are known to ex- 
hibit, for nearby trajectories, a subexponcntial deviation 
of the type 5xt ~ 8xg exp(A Q i Q ), where a = l/(z — 1). 
According to the Aaronson-Darling-Kac (ADK) theorem 
[2(i| , for randomly distributed initial conditions and suffi- 
ciently large times, the ratio X a / (A), where (A) is a suit- 
able average for the exponents A a , converges in distribu- 
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FIG. 5: Curve (a): probability of having X a < (A) in Pomeau- 
Manneville maps as function of a £ (0, 1), according to Eq. 
(|26p . The marked point (a ~ 0.747) corresponds to the 
equiprobability. Detail: curves (b) and (c) are, respectively, 
the probability density of Mittag-Leffler distribution with unit 
first moment and a — 0.1 (left j/-scale) and a = 0.9 (right y- 
scale) in the interval < x < 3. It is clear why the probability 
decays for increasing a: the probability density tends towards 
a 5-function centered in x — 1, but from the right handed side. 
(See curve (c) in Fig. [3]). 



tion terms towards a Mittag-Leffler random variable with 
unit first moment and index a £ (0,1). Such statistics 
was also considered previously in Ref. [18j from the nu- 
merical point of view. Some recent numerical works (l9j 
have reported a regular tendency of A Q be smaller than 
the average (A) for large values of z (small a). In fact, 
in the first moment (A) is calculated differently from 
the ADK theorem, it is obtained there from a continuous- 
time stochastic model, but its values are, for the con- 
sidered Pomeau-Manneville maps, the same of the ADK 
ones. Since X a / (A) is a random Mittag-Leffler variable, 
we can evaluate the probability of having A Q < (A) 

Prob(A Q < (A)) = f p { :\x)dx = 1 - A Q (r 1 / Q (l + a)), 
Jo 

(26) 

where r = r 1 / Q (l + a) assures that pa\x) has unit fist 
moment, as required by the ADK theorem. Fig. [5] de- 
picts this probability as function of a. For a — > 0, we 
have Prob(A Q < (A)) -> 1 - 1/e « 63%. As we see, 
the tendency reported in [l9[ of having A Q < (A) can be 
clearly understood from the ADK theorem. Moreover, 
the aspect of the Mittag-Leffler distributions for small a 
and a — > 1 explains why these intermittent systems do 
exhibit such kind of behavior. For a — > 0, the probability 
density function of the Mittag-Leffler distribution has the 
form p a (x) s» e~ x (See Fig. [5]). Its maximum is located 
at x — 0, and it is clear that the typical values of the ran- 
dom variable are always smaller than its average. On the 
other hand, for a — > 1 (z — > 2 + ), the probability density 
resembles a S- function with center approaching x = 1, 
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but from the right handed sided, see Fig. In this case, 
the typical values of the random variable remain close the 
value of its average. Also from figure, we have that for 
a > 3/4 the probability of having A Q > (A) is favorable 
over Eq. (f26|) . (In fact, the equiprobability corresponds 
to aw 0.747.) In terms of the distribution of Lyapunov 
exponents for the Pomeau-Manneville maps, this would 
correspond to have a slight predominance of A Q greater 
than the average (A) for small z > 2. This seems, in fact, 
marginally evident from (ljjj . but further work is neces- 
sary to establish this fact with the same certainty of the 
behavior for large z. This kind of problem in intermit- 
tent systems are very interesting and certainly deserve a 
deeper investigation. 

IV. FINAL REMARKS 

Motivated by the closed form formulas in terms of the 
Meijer G-functions for the probability densities g a {x) of 
one-sided Levy distributions with rational a — l/k pro- 
posed by Penson and Gorska in |3j, we have introduced 
a numerical scheme for the computation of g a (x) for any 
real a £ (0,1). By exploring the relation between one- 
sided Levy and Mittag-Leffler distributions, we extend 



our procedures to include the evaluation of the the proba- 
bility densities p a (x) of Mittag-Leffler distributions. The 
main advantage of our numerical approximation is that it 
can be applied for any value a, while Penson and Gorska 
formula ([2]) is rather problematic for fractions with large 
denominators. As an application of our procedures, we 
determine the maximum location and value for the den- 
sities g a (x) and p a {x) as function of the index a € (0, 1). 
We show that a « 0.567 and a « 0.605 correspond, re- 
spectively, to the one-sided Levy and Mittag-Leffler dis- 
tributions with shortest maxima. Furthermore, we use 
our numerical procedure for the evaluation of Mittag- 
Leffler distribution to show that a recently described sta- 
tistical behavior for intermittent systems [19| , namely the 
predominance of having Lyapunov exponents A a smaller 
than the theoretical average (A) for Pomeau-Manneville 
maps with large z, is nothing else than a consequence 
of the Mittag-Leffler statistics. We hope our numerical 
procedures could be useful for this kind of study. 
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